#include <iostream>
#include <fstream>
#include <vector>
#include <string>

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

#include <iostream>
#include <fstream>
#include <typeinfo>

#include "util.h"
#include "diffuse.h"

// 方阵的维度，
#define nMAX 3

using namespace std;

int main()
{

	char str1[64], str2[64], str3[64];

	sprintf(str1, "output_data");
	if (!access(str1, 0) == 0)
		mkdir(str1, 0777); // if  no "data" folder, construct it for output date

	char buf_data[40];

	// double  a[nMAX][nMAX] = { 0 };//定义一个nMAX*nMAX*的矩阵，用于存放数据
	double **a = dmatrix(1, nMAX, 1, nMAX); // 或者我们可以动态开辟数组，这种方式记得在程序结束后释放内存。使用方式 free_dmatrix(a, 0, nMAX, 0, nMAX);
	double *b = dvector(1, nMAX);

	// 怎样读取数据文件，比如data.txt中存储的是矩阵A数值
	ifstream infile, infile2;		// 定义读取文件流，相对于程序来说是in
	infile.open("data.txt");		// 打开文件
	for (int i = 1; i <= nMAX; i++) // 定义行循环
	{
		for (int j = 1; j <= nMAX; j++) // 定义列循环
		{
			infile >> a[i][j]; // 读取一个值（空格、制表符、换行隔开）就写入到矩阵中，行列不断循环进行
		}
	}
	infile.close(); // 读取完成之后关闭文件

	infile2.open("data_b.txt");		// 打开文件
	for (int i = 1; i <= nMAX; i++) // 定义行循环
	{
		infile2 >> b[i]; // 读取一个值（空格、制表符、换行隔开）就写入到矩阵中，行列不断循环进行
	}
	infile2.close(); // 读取完成之后关闭文件

	// LU 分解函数
	double **L = dmatrix(1, nMAX, 1, nMAX); // 开辟内存，用存储LU分解后的下三角形矩阵L
	double **U = dmatrix(1, nMAX, 1, nMAX); // 开辟内存，用存储LU分解后的上三角形矩阵U
	double *x = dvector(1, nMAX);

	//LU_decompose(a, L, U); // 算法2.2
	Gauss_solver(a, x, b, nMAX); // 算法2.4 double **a, double *x, double *b, int n

	sprintf(buf_data, "output_data");
	//print_data(a, L, U, buf_data);
	print_data2(a, x, buf_data);

	// 动态开辟的内存，程序结束前，最好释放掉， 不做内存释放，程序运行后，系统也会自动释放内存
	free_dvector(x, 1, nMAX);
	free_dvector(b, 1, nMAX);
	free_dmatrix(a, 1, nMAX, 1, nMAX);
	free_dmatrix(L, 1, nMAX, 1, nMAX);
	free_dmatrix(U, 1, nMAX, 1, nMAX);
	return 0;
}

/* 输入变量：a
输出变量L， U*/
// 算法2.2
void LU_decompose(double **a, double **L, double **U) // 讲义算法2.2
{
	zero_matrix(L, 1, nMAX, 1, nMAX); // 矩阵使用前，数据清零
	zero_matrix(U, 1, nMAX, 1, nMAX); // 矩阵使用前，数据清零

	// 将L设置为单位矩阵
	for (int i = 1; i <= nMAX; i++)
	{
		for (int j = 1; j <= nMAX; j++)
		{
			if (i == j)
			{

				L[i][j] = 1;
			}
			else
			{
				L[i][j] = 0;
			}
		}
	}

	int i, j, k;

	for (k = 1; k <= nMAX-1; k++)
	{
		for (j = k; j <= nMAX; j++)
		{
			// cout<<k<<j<<endl;
			U[k][j] = a[k][j];
		}
		for (i = k + 1; i <= nMAX; i++)
		{
			L[i][k] = a[i][k] / a[k][k];
			for (j = k + 1; j <= nMAX; j++)
			{
				a[i][j] = a[i][j] - L[i][k] * U[k][j];
			}
		}
	}
	U[nMAX][nMAX]=a[nMAX][nMAX];

	return;
}

void LU_row(double **a) // 讲义算法2.4  行优先存储，L和U直接存储在原来的a矩阵中
{
	int i, j, k;
	for (i = 2; i <= nMAX; i++)
	{
		for (k = 1; k <= i - 1; k++)
		{
			a[i][k] = a[i][k] / a[k][k];
			for (j = k + 1; j <= nMAX; j++)
			{
				// cout << k << i << j << endl;
				a[i][j] = a[i][j] - a[i][k] * a[k][j];
			}
		}
	}

	return;
}

// Ax=b, dim: n
void Gauss_solver(double **a, double *x, double *b, int n) // 讲义算法2.4  行优先存储，L和U直接存储在原来的a矩阵中
{
	int i = 0, j = 0;

	double *y = dvector(1, n);
	for (int i = 1; i <= nMAX; i++)
	{
		y[i] = 0.0;
	}

	LU_row(a);
	double **L = dmatrix(1, nMAX, 1, nMAX);
	for (int i = 1; i <= nMAX; i++)
	{
		for (int j = 1; j <= nMAX; j++)
		{
			if (i > j)
			{
				L[i][j] = a[i][j];
			}
			else if (i == j)
			{
				L[i][j] = 1;
			}
			else
			{
				L[i][j] = 0;
			}
		}
	}

	y[1] = b[1] / L[1][1];
	for (i = 2; i <= n; i++)
	{
		for (j = 1; j <= i - 1; j++)
		{
			b[i] = b[i] - L[i][j] * y[j];
		}
		y[i] = b[i] / L[i][i];
	}

	x[n] = y[n] / a[n][n];

	for (int i = n - 1; i >= 1; i--)
	{
		for (int j = n; j >= i + 1; j--)
		{
			y[i] = y[i] - a[i][j] * x[j];
		}
		x[i] = y[i] / a[i][i];
	}
	for (int i = 1; i <= n; i++)
	{
		cout << x[i] << endl;
	}

	free_dvector(y, 1, n);

	return;
}

void Gauss_partial_pivoting(double **a, double *x, double *b, int n) // 讲义算法2.8 // 部分选主元高斯消去法
{
	int i, j, k, q;

	int *p = ivector(1, n);

	for (i = 1; i <= n; i++)
		p[i] = n; // 生成1到n的自然数列

	double amax;
	double atem;
	int ptem;

	for (k = 1; k <= n - 1; k++)
	{
		for (i = k; i <= n; i++) // 选列主元，其中q表示主元所在的行
		{
			amax = a[i][k];
			q = i;
			if (amax >= a[i][k])
			{
				amax = a[i][k];
				q = i;
			}
		}

		if (q != k) // 交换a的第k行与第q行
		{
			atem = a[k][j];
			a[k][j] = a[q][j];
			a[q][j] = atem;
		}

		ptem = p[k];
		p[k] = p[q];
		p[q] = ptem; // 更新置换矩阵

		for (i = k + 1; i <= n; i++)
		{
			a[i][k] = a[i][k] / a[k][k];
			for (j = k + 1; j <= n; j++)
			{
				a[i][j] = a[i][j] - a[i][k] * a[k][j];
			}
		}
	}

	free_ivector(p, 1, n);

	return;
}

void print_data(double **phi, double **eta, double **xx, char *buf)
{
	extern long nx, ny;

	char buffer_phi[32], buffer_eta[32], buffer_xx[32];

	FILE *fphi;
	FILE *feta;
	FILE *fxx;

	sprintf(buffer_phi, "./%s/a.m", buf);
	sprintf(buffer_eta, "./%s/L.m", buf);
	sprintf(buffer_xx, "./%s/U.m", buf);

	fphi = fopen(buffer_phi, "w");
	for (int i = 1; i <= nMAX; i++)
	{
		for (int j = 1; j <= nMAX; j++)
		{
			fprintf(fphi, " %16.14f ", phi[i][j]);
		}
		fprintf(fphi, "\n");
	}
	fclose(fphi);

	feta = fopen(buffer_eta, "w");
	for (int i = 1; i <= nMAX; i++)
	{
		for (int j = 1; j <= nMAX; j++)
		{
			fprintf(feta, " %16.14f ", eta[i][j]);
		}
		fprintf(feta, "\n");
	}
	fclose(feta);

	fxx = fopen(buffer_xx, "w");
	for (int i = 1; i <= nMAX; i++)
	{
		for (int j = 1; j <= nMAX; j++)
		{
			fprintf(fxx, " %16.14f ", xx[i][j]);
		}
		fprintf(fxx, "\n");
	}
	fclose(fxx);
	cout << "output ok" << endl;
	return;
}
void print_data2(double **phi, double *x, char *buf)
{
	extern long nx, ny;

	char buffer_phi[32], buffer_eta[32];

	FILE *fphi;
	FILE *feta;

	sprintf(buffer_phi, "./%s/a.m", buf);
	sprintf(buffer_eta, "./%s/x.m", buf);

	fphi = fopen(buffer_phi, "w");
	for (int i = 1; i <= nMAX; i++)
	{
		for (int j = 1; j <= nMAX; j++)
		{
			fprintf(fphi, " %16.14f ", phi[i][j]);
		}
		fprintf(fphi, "\n");
	}
	fclose(fphi);

	feta = fopen(buffer_eta, "w");
	for (int i = 1; i <= nMAX; i++)
	{
		fprintf(feta, " %16.14f ", x[i]);
		fprintf(feta, "\n");
	}
	fclose(feta);

	cout << "output ok" << endl;
	return;
}
